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We demonstrate that Widom's particle insertion technique provides a convenient and efficient method to determine the effective 
pair interaction between complex, composite soft-matter particles in the zero-density limit. By means of three different test 
systems, i.e. amphiphilic dendrimers, electrostatic polymers and colloids coated with electrostatic polymers, we demonstrate the 
validity and the power of the presented method. 



1 Introduction 

Soft materials, be they colloids, polymers or proteins, are of- 
ten complex constructs that live in a bath of small molecules. 
In addition, these mesoscopic particles themselves are often 
composite objects that contain flexible moieties^] As a result, 
atomistic modeling of the structure and dynamics of soft ma- 
terials would need to span a wide range of length and time 
scales, which makes such an approach infeasible in all but the 
simplest cases (few mesoscopic particles and short times). 

In an effort to simplify the description of the system at hand, 
coarse-grained models are being developed that aim to capture 
the mesoscopic and macroscopic behavior of the system by 
including in the description only the most relevant degrees of 
freedom. In the simplest coarse-graining approach each com- 
plex, mesoscopic particle is characterized by only one effec- 
tive coordinate. In the case of particles that have on average 
inversion symmetry, such as e.g. polymer-coated colloids, this 
is the coordinate of the center of inversion symmetry. For sys- 
tems without such symmetry (or, to be more precise, where 
this symmetry center cannot be identified with a fixed coor- 
dinate in the molecular frame) it is conventional to use the 
center of mass to specify the position of the particle, which is 
the procedure that is followed for polymers or proteins. The 
mesoscopic particle is then represented as a "soft" sphere and 
the effective interaction potential determines the softness. The 
simplest effective potentials are obtained in the low-density 
limit where many-body interactions can be ignored. In spite 
of their simplicity, such models have been shown to work suc- 
cessfully for a range of mesoscopic systems (see e.g.W^j). 

Computing effective interactions for arbitrarily complex 



t In what follows, we use the term "mesoscopic" to describe particles that con- 
tain internal degrees of freedom that will be integrated out in a more coarse- 
grained description. Usually, these particles will be in the mesoscopic size 
range from nano- to micrometers. 
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mesoscopic particles is one of the key steps in the develop- 
ment of a coarse-grained model. There are many ways in 
which these interactions can be obtained. The present paper 
describes an approach that we found to be considerably more 
efficient than other approaches that we explored. 

The paper is organized as follows: in Sec. [2] we define the 
effective interaction and explain standard Monte Carlo (MC) 
techniques to determine it, highlighting their strengths and 
their limitations. In Sec.|3]we adapt Widom's particle insertion 
method to the present problem and show that it is an efficient 
way to determine effective interactions. In Sec. [4] we explain 
how to implement this formalism and test the method for three 
different model systems in Sec. [5] comparing the results from 
the different techniques. In the concluding section (Sec.|6| we 
point out possible limitations of the method that we present. 

2 Computing effective interactions: the stan- 
dard approach 

Consider two mesoscopic particles confined in a volume V and 
positioned at a distance R\2 between them. At infinite dilution, 
the effective pair interaction ^ e fi(Rn) between these particles 
is related to their radial distribution function^ g(R 12) via 



g(*i2) = eatp[-P* eff (^ u )]. 



(1) 



where p = 1/kT is the reciprocal temperature 59 . The radial 
distribution function, g(Rn), can easily be measured during 
the simulation and has to be normalised to 1 for large separa- 
tions between the particles. Eqn. [T]thus offers a straightfor- 
ward way to measure the effective interaction within simula- 
tions. However, the repulsion between two mesoscopic par- 
ticles usually increases as the particles approach each other, 
possibly reaching several kT and configurations where parti- 
cles are that close are rare. As a result, the relative error in 
g(Rn) at short distances will therefore be large, which results 
in a concomitantly large error in < t > e ff(^i2)- In order to ob- 
tain an accurate estimate of 't> e g(Ri2) without wasting time on 



irrelevant, though easily accessible configurations, it is there- 
fore necessary to use a simulation technique that samples also 
the regime where <Peff(Rn) 3> kT. 

One way to overcome this sampling bottleneck is to sample 
all relevant values of R12 more or less uniformly. Conceptu- 
ally, the simplest approach to achieve this is to introduce an 
external biasing potential that acts on the effective coordinates 
and counteracts the repulsion between the particles. The total 
interaction is then given by <P(Rn) — ^effC^n) + ^1^(^12) 
and the radial distribution function is modified to read 



g(R n ) - exp[-p*(/?i 2 )] 

= exp[-P4> eff (^i 2 )]exp[ 



P*bias(^12) 



Comparing this last equation to Eqn. [T] we see that to deter- 
mine the radial distribution function of the unbiased system, 
g(Ri2)> we have to increment the histogram of the probabil- 
ity of finding the molecules separated by a certain distance 
R12 in each measurement by exp [p<J>bi as (^i2)] rather than by 
1 . The effective interaction can then be determined from this 
histogram via Eqn. [T] 

However, there are two major disadvantages to this ap- 
proach. First, to sample all relevant distances uniformly, 
<P(Rn) should be zero (or constant) for all R\2 and there- 
fore the biasing potential should ideally be ObiasC^^) = 
— <P e ff(Ri2), which would correspond to knowing a priori the 
sought-after answer. In practice, long simulations are often 
needed where the biasing potential is determined in a trial- 
and-error procedure, enhancing the guess for ObiasC^n) in 
a rather cumbersome iterative process until the difference to 
— <J> e ff(^i2) is smaller than kT, which will then allow for an 
efficient sampling of all distances. But even once the bias- 
ing potential is known with sufficient accuracy, sampling will 
still be slow due to the fact that the mesoscopic particles have 
to diffuse over the relevant range of R12 values to visit each 
distance often enough. 

In an automated process, one can use the Wang-Landau ap- 
proach 10 to determine the biasing potential and thereby the ef- 
fective interaction in an iterative way. In this scheme, g{Rn) is 
initially unknown and set to unity for all distances R\2- During 
the simulation, moves from an old distance R° 2 with potential 
energy U(R" 2 ) to a new distance R" 2 with potential energy 
U(R" 7 ) are accepted with a modified Metropolis acceptance 
rule, where the probability to accept the move is given by 



P lux (R u ^R n 12 )=rmn 



In each measurement, the histogram for g at the current posi- 
tion R12 is multiplied by a factor /, which is usually taken to 
be 2. At the same time, a histogram of the distances visited is 
measured and once this histogram is sufficiently flat, a coarse 
guess for g(Rn) has been obtained. To refine the results, / is 



set to \ff, and the simulation is iterated like this until / ~ 1. 
Then, the effective interaction can be determined via Eqn.[T] 

Alternatively, one can break down the range of distances 
into several small windows^ 11 ! 12 ! and then use simple bias- 
ing potentials to force the system to stay within each win- 
dow. Typical choices for these so-called umbrella potentials 
are e.g. hard walls^ or spring-like potentials "t^iasC^ 12 ) = 
\kj(R\2 — Rj) 2 , where the kj are spring constants that deter- 
mine the width of each window j located at Rj. Carrying out 
separate simulations for the different windows, one can sys- 
tematically vary the separations between the mesoscopic par- 
ticles. Within each of these windows j, a separate histogram 
of the probability of finding the molecules a certain distance 
apart, g^(R), is recorded. At the end of the simulations, the ef- 
fective potential Oeff (^12) between the molecules is obtained 
by merging the effective potentials obtained within each win- 
dow, 



12 



■ln[gJ(R l2 )]-^L(Rn) + Cj 



where Cj is a normalization constant. Since the various Cj 
are initially unknown, the concatenation of < & e ff(Ri2) will dis- 
play discontinuities at the windows' edges. To obtain a con- 
tinuous 4 ) e ff(/?i2), the cj are chosen such that the data are 
aligned to each other at the edges of the windows or, more 
sophist icatedly , the multiple histogram method can be imple- 
mented ^ 14 l 15 [ Using this umbrella sampling, care has to be 
taken to choose the windows and umbrella potentials in a way 
that the variation of the effective interaction within each win- 
dow does not exceed 1-2 kT to reliably sample the entire win- 
dow. Moreover, to improve quality of the matching of the 
different parts of < t > e ff(^i2) at the edges of the windows, it is 
advisable to use overlapping windows. 

In view of the various limitations of above techniques, it 
would be desirable to have a straightforward, unbiased method 
of sampling all distances of approach between two mesoscopic 
particles with the same probability, while obtaining an effi- 
cient estimate of the effective interaction. We show in the next 
section how Widom's particle insertion method can be used to 
achieve precisely this. 

3 Formalism 

The following formalism is based on Widom's insertion 
method to calculate the chemical potential of non-uniform 
fluids^. We adapt his approach to determine the effective 
interactions between two mesoscopic particles in the zero 
density limit in a fast and straightforward manner. The 
method presented here is also a generalisation of the scheme 
used in Ref. 2, where the effective interaction between two 
self-avoiding, i.e. hard-core, polymers was determined from 
the overlap probability when placing the polymers random 
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distances apart. Despite this generalisation being straight- 
forward, we are not aware of earlier implementations. 

We consider a fluid of N composite, mesoscopic particles 
confined in a volume V. The effective coordinate of particle i 
is denoted by vector Rj with i = 1 , . . . , N. 

The probability P to find this system in a certain configura- 
tion {Ri , . . . ,Rjv} is given by the Boltzmann distribution 



P({R u ...,R N }) = —e 



-P£/(R, 



U is the potential energy of the system including all intra and 
inter-particle interactions. Zn is the configurational integral 
given by 



dRL.-dRve-™ 



As Widom argued^, the density profile at position r within 
the volume can be written as 



P(r) 



N 



dR, ...dR 



-pf/(Ri,...,R w _i,r) 



Splitting the potential energy in above equation into a contri- 
bution due to all the interactions in the (N — 1) -particle system 
and an energy change, At/, due to adding an Mh particle at 
position r, i.e. 

l/(Ri,...,R A ,_i,r) = l/(Ri,...,R w _i)+AU(Ri ) ...,H JV _i,r) 
we can write that 

p(r) = ^ f dR!...dRiv-ie- p£/(Rl «w-i) 4f -PAi7(Ri,..-^-i^) j 
Zn J 

which simplifies to 



P(r) 



= NZ N -i i _ pAt/(Rl Rn l T 



N-l.r 



where (x) 



N-l,r 



denotes an ensemble average over quantity x 



in the (N — 1) particle system at position r. 

Now, we apply this last formula to a system of only two 
mesoscopic particles, labeled 1 and 2, whose effective coordi- 
nates are positioned at Ri and r = R2. Then, 



P(R 2 ) 



NZi 



-PAU(Ri,R 2 ) 



N=1,R 2 



Since we want to determine a radially symmetric interaction 
potential, we can position the first particle in the origin and 
place the second at a distance R12 = |R2 — Ri | = |R2 1 - The 
above equation then simplifies to 



P(*12) 



NZi 



The density profile can be expressed as p(Rn) = Pg(Rn), 
where p is the average number density. Recalling Eqn. [T] we 
find that 



P<£effORi 2 ) = -ln 



,-$AU(R l: 



N=1,R [2 



(2) 



The constant c cannot be directly measured since it is propor- 
tional to the configurational integral of the system. However, it 
is expected that the range of the effective interaction between 
two typical macromolecules is finite, i.e. <t> e ff — > as R — s- °°. 
Therefore, with /? max being a separation between the meso- 
scopic particles where their interaction has decayed to zero, c 
can be determined as 
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4 Simulation Technique 

Eqn. [2] allows us to determine the effective interaction be- 
tween two mesoscopic particles as a function of distance R\2 
between their effective coordinates R;, i = 1,2, in a straight- 
forward, efficient and unbiased way. 

Let the coordinates of the M constituent parts of the meso- 
scopic particles be given by {r,- 1, . . . ,r,jw}, i — 1,2. In a 
first step, we equilibrate each of the two mesoscopic parti- 
cles in isolation and sample their configurational space with 
a standard Metropolis MC simulation. The particle moves 
employed will depend on the particular system under study 
and will be chosen to ensure ergodicity. Each MC move 
from an old internal configuration o, = {r"j , . . . , r" M } to a new 
one rit — {if 1, ■ ■ ■ i = 1(2, is accepted according to the 

Metropolis MC acceptance rule 17 , i.e. with a probability 



^acc [pi 



min ( 1 , e 



N=l,R l2 



where U (x,) is the intramolecular potential energy of configu- 
ration x of particle i. 

Once the mesoscopic particles have been equilibrated, we 
randomly sample from the equilibrated conformations that 
are used in the second step of the algorithm. We fix the 
coarse-grained coordinate of the first particle at the origin. 
We then generate effective coordinates of the second parti- 
cle, uniformly distributed in the interval R\2 G [0,R max } and 
we measure the intermolecular potential energy AU (Rn) be- 
tween the two mesoscopic particle^] This procedure yields 
e -$&u(Rn) as a function of R\2 with high efficiency, since no 
acceptance test has to be applied in this last step. Repeated 

% To ensure good statistics at close approach of the mesoscopic particles, we 
sample all distances uniformly on a line segment. An alternative, but less 
accurate procedure would sample points uniformly distributed in a spherical 
volume, but such a procedure would disfavor the sampling of short distances. 
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sampling of steps one and two allows us to determine the av- 
erage in Eqn. [2] and thereby the effective interaction between 
the two molecules. 

5 Application to test systems 

To demonstrate the power and validity of the method pre- 
sented above, we determine the effective interactions of three 
different test systems. 

5.1 Electrostatic polymers 

As a first, simple test system, we study a strongly coarse- 
grained representation of single-stranded DNA 18 , where the 
strands are represented as electrostatically charged, freely 
jointed chains. Vertices are separated by segments of fixed- 
length /Kuhn = l-5nm 19 and interact with each other via a 
Debye-Hiickel potential which reflects both the charge of the 
sugar-phosphate backbone of the DNA strands and the solvent 
conditions 18 , 



2 exp(-Kr , '>) 

: i b n — - 



where r' J is the distance between vertices i and j, K is the in- 
verse Debye screening length, q corresponds to the charge per 
vertex and lg is the Bjerrum length. Here, we choose k _1 = 
0.67nm, lg = 0.7nm and q = 3.05e^and sample chains of 20, 
30, 40 and 50/Kuhn with both biased simulations and Widom's 
particle insertion technique. In the former, we sample the 
polymers for 1 x 10 7 MC sweeps, where each sweep consists 
of 500 MC moves (i.e. either a pivot mov d 21 ! 2 ^ , crankshaft 
moveB^l rotation or translation of a polymer). For Widom's 
method, we decorrelate each polymer for 500 MC moves, us- 
ing pivot and crankshaft moves, and with the two configura- 
tions thus obtained we sample 500 different distances between 
the polymers. We then repeat this cycle 2.5 x 10 6 times. In 
Fig. [T] we compare the results from the two methods and 
see that they perfectly coincide, confirming the validity of 
Widom's method. 

To show how rapidly Widom's method converges, we com- 
pare to the results from a brute force simulation and a biased 
one with Obias = — ^eff- We ran all three simulations for the 
same amount of computational time, i.e. 1 hour on a proces- 
sor of our quad-core dual Xeon (Harpertown) cluster. For the 
Widom method, this corresponds to 6 x 10 4 cycles, for the 
brute force and the biased simulation to 1 1 x 10 4 MC sweeps 
each. As shown in Pig. [2j the results from Widom's insertion 
are equilibrated already after this little of simulation time. By 
contrast, the brute force simulation has not yet visited small 
distances at all and the data of the biased simulation shows 
considerable noise. This impressively highlights the efficiency 
of Widom's method to determine effective interactions. 




Fig. 1 (Color online) Comparison between Widom's method (solid 
lines) and biased simulations (symbols) of the effective interaction 
3> e ff between two electrostatically charged polymer chains of length 
20 (diamonds), 30 (circles), 40 (squares) and 50/Kuhn (red, triangles) 
as a function of the distance between their centers of mass. The inset 
shows a simulation snapshot of two chains of 20/Kuhn- 



5.2 Colloids coated with electrostatic polymers 

As a second example, we study a system of colloids coated 
with electrostatically charged polymers that can move on the 
colloid's surface. In this case, we determine the effective in- 
teraction between the centers of the colloids rather than the 
centers of mass. The vertices of the chains interact amongst 
each other as described in Sec. |5.1| while the colloids are im- 
penetrable spheres of fixed radius R co \\. Again, we implement 
crankshaft and pivot moves for the polymer chains, but we 
also regrow them completely via a configurational bias MC 
algorithm^, using 5 trial directions in each step of the re- 
growth. The colloids are subject to rotations and translations. 
Our model system is made of colloids of Z? co rj = 4/Kuhn> each 
coated with 10 chains of 5 Kuhn segments. For Widom's 
method, we simulated 4 x 10 6 cycles consisting of 250 MC 
steps to decorrelate the coating of the colloids before sampling 
500 different distances between the colloids. 

In Fig. [3] we see that the Widom particle insertion data per- 
fectly reproduces the results of unbiased, brute force simula- 
tions, which we ran for 1 x 10 7 MC sweeps. In the latter, each 
sweep consisted of 220 MC moves (i.e. pivot or crankshaft 
move, configurational bias MC regrowth, colloid rotation or 
translation). While the brute force simulations needed two 
days of simulation time and were averaged over 6 different 
simulation runs to obtain reliable statistics, the data from a 
single simulation applying Widom's method already shows 
decent statistics after as little as 5 x 10 4 cycles, which cor- 
responds to 1 .5 hours of simulation time on a single processor 
of our cluster. 
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Fig. 2 (Color online) Comparison between Widom's method (solid 
line), biased simulations (solid line with circles) and brute force 
calculations (solid line with squares) of the effective interaction <J> e ff 
between two electrostatically charged polymers of length 20 after 
one hour of simulation time. While the data from Widom's insertion 
is equilibrated, statistics on the data of the other two methods is still 
insufficient. 



Fig. 3 (Color online) Comparison between Widom's method (solid 
line) and brute force simulations (crosses) of the effective 
interaction <I> e ff between two colloids of radius R co \\ = 4/Kuhn coated 
with 10 electrostatically charged polymer chains of length 5/Kuhn as 
a function of the distance between their centers. The inset shows a 
simulation snapshot of the colloids. 



5.3 Amphiphilic dendrimers 

Our final example concerns the determination of the effective 
interaction between two amphiphilic dendrimers 624 , which 
are expe cted to show clustering behavior at sufficiently high 
densities02H3. 

We sample two different second-generation amphiphilic 
dendrimers. They have two central monomers and while the 
end-groups form a solvophilic shell, all inner monomers repre- 
sent the solvophobic core. The bonds between monomers are 
modeled by the finitely extensible nonlinear elastic potential, 
while all other intera ctions b etween monomers are modeled 
by the Morse potential ^ 28 ! 29 !. For the parameters of the differ- 
ent interactions, we choose the same as for dendrimer D5 and 
D-i from Ref. 30 and we refer to this reference for details. 

Again, we implement Widom's insertion method and com- 
pare our results to the effective interactions found in Ref. [30] 
via umbrella sampling. For the umbrella sampling, 15 slightly 
overlapping windows with /? max = 5R g were used, where R g 
is the radius of gyration of a single dendrimer. The systems 
were then sampled for 2 x 10 8 MC sweeps in each window. 
These simulations can be carried out in parallel. For Widom's 
method, we sampled the dendrimers for 4 x 10 6 MC cycles. 
In each of these cycles, each dendrimer was decorrelated in 
280 single monomer displacements and the resulting configu- 
rations were used to sample 1000 different distances between 
them. 

As can be seen from Fig. [4] Widom's method manages to 
capture the interaction potential even at close approach of the 
dendrimers where deformations to the dendrimer' s conforma- 



tions are expected. While Widom's technique already gives 
reliable statistics for the effective interaction after as little as 

5 x 10 5 cycles, corresponding to 3h of simulation time, um- 
brella sampling had to be carried out on 15 processors for 
roughly a day to collect the necessary statistics. 

6 Conclusions 

In this contribution we have shown that Widom's particle in- 
sertion method can be adapted to determine the effective in- 
teraction between two mesoscopic particles in a completely 
unbiased, efficient way. While we showed the success of the 
method for three rather different model systems, this method 
can not be applied in its present form to particle insertion in 
an explicit solvent. Further, we expect the method to fail for 
those systems where close approach between the two particles 
leads to a very strong deformation of the molecules' confor- 
mations that will not be sampled adequately in the ideal-gas 
reference system that we use to generate independent confor- 
mations. The usual way to test for the reliability or the failing 
of Widom's method is the so-called overlapping distribution 
method^ED where two simulations are needed: one of a one- 
particle system, inserting a second one at frequent intervals; 
and one of a two-particle system, removing one at frequent in- 
tervals. Widom's method will only give reliable answers when 
the probability distributions of finding certain energy differ- 
ences AU upon insertion and removal, respectively, have suf- 
ficient overlap^. However, in these cases Widom's method 
will at least allow for an educated guess of a biasing potential 
that can be used with one of the standard techniques, and will 
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Fig. 4 (Color online) Comparison between Widom's method (solid 
lines) and brute force simulations (symbols) of the effective 
interaction 4> e ff between two amphiphilic D5 (circles) and D7 
(squares) dendrimers as a function of the distance between their 
centers of mass. The inset shows a simulation snapshot of two 
dendrimers. 



thereby allow for a considerable speed-up of the conventional 
methods. 

In summary, the application of Widom's particle insertion 
method to the determination of effective interactions between 
two mesoscopic particles has several advantages over the stan- 
dard techniques: first, it allows for unbiased simulations, 
while the standard techniques require a guess of appropriate 
biasing potentials to force the system into close approaches. 
These potentials have to be determined with high accuracy in 
an iterative process, or if less precise potentials are used, simu- 
lations in multiple windows of different separation ranges are 
required. Furthermore, Widom's method is very easy to im- 
plement. Finally, biased or brute force simulations rely on the 
diffusion of the particles through the simulation box, where 
each change in separation will only be accepted with a cer- 
tain probability depending on the ratio between the Boltzmann 
weights of the old and the new configurations. This makes the 
relative diffusion of highly entangled particles slow and can 
lead to slow statistical fluctuations in the data of the effective 
interaction. By contrast, within Widom's method particles are 
simply placed at a given separation and the intermolecular en- 
ergy difference is measured, without having to employ any 
acceptance/rejection step. Therefore, Widom's particle inser- 
tion method immediately allows for an exploration of even 
the closest approach and will give reliable statistics for the 
effective interaction for all ranges of separation within very 
short simulation times. For the cases that we have studied, 
this makes Widom's particle insertion method the better tech- 
nique to compute effective interactions between mesoscopic 
particles. 
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